
% Read both images
M=im2double(imread('lenag1.png'));
F=im2double(imread('lenag2.png'));
% Parameters
options.sigma_fluid=8; options.sigma_diff=0; options.alpha=1.5;
% Transformation field
T=zeros([size(M) 2]);
% Store the dimensions of transformation field, and make, a long vector from T
sizes=size(T); T=T(:);
% Start the demon energy registration optimizer
optim=struct('Display','iter','GradObj','on','MaxIter',500,'OutputFcn', @store_transf);
T=fminlbfgs(@(x)demons_energy(M,F,[],[],x,sizes,options),T,optim);
% Reshape O_trans from a vector to a matrix
T=reshape(T,sizes);
% If diffusion regularization is used: smooth transformation field
if(options.sigma_diff>0)
    Hsmooth=fspecial('gaussian',[options.sigma_diff*6 options.sigma_diff*6],options.sigma_diff);
    T(:,:,1)=imfilter(T(:,:,1),Hsmooth); T(:,:,2)=imfilter(T(:,:,2),Hsmooth);
end
% Transform the input image
Ms=movepixels(M,T(:,:,1),T(:,:,2));
% Show the registration results
figure,
subplot(1,3,1), imshow(M); title('input image 1');
subplot(1,3,2), imshow(F); title('input image 2');
subplot(1,3,3), imshow(Ms); title('transformed image 1');